%% An example for computing the EM flux

clear all

I = 1;
dw = 0.32004e-3; %%28AWG wire
Area    = dw^2;
J     = I/Area; 
nu0 = 4*pi*1e-7;

% % %%EM4
a_o    = 15.875e-3;
a_in   = 4.76e-3;
l     = 9.525e-3;

em_n = 16;
em_k = 6;
em_delta_l = 0.442;
% error = 0.0428;
em_l_bar = l*em_delta_l;
m_em =   1.0e-009 *[120.3 44.6 131.89 133.96 134.85 108.01 48.25];

%% find the equivalent single layer model for the EM
ro = 2*a_o/l;
ri = 2*a_in/l;

Z = 0.5;
rou = fzero(@(r)s_radius(ro,ri,r,Z),0.5);
ae = l*rou/2;
B3d = J*MagneticField3D(0,0,l/2,l,a_in,a_o);   %%coil
B_esl = MagneticField_ESL(0,0,l/2,l,1,ae);
Je_dw = B3d(3)/B_esl(3)

%% compute the flux useing ESL, ML, and DMP for the EM
[m_Nem,p_emt,p_emb] = pole_cyl(m_em,a_o,em_l_bar,0,em_n,em_k);
p = [p_emt,p_emb];
m = [m_Nem,-m_Nem];

%%along radius
% r = 0:2e-2:2*a_o;
r = 0:1e-4:l;
n = size(r,2);
np = size(p,2);
e = 0.2e-3;

% %%along z axis
% r = l/2:2e-5:1.5*l;
% n = size(r,2);
% np = size(p,2);

for j = 1:n
%     along radius
%     x_d = r(j);
%     y_d = 0;
%     z_d = l/2+e;
%     along z axis
    x_d = 0;
    y_d = 0;
    z_d = r(j);
    
    B_ml = MagneticField3D(x_d,y_d,z_d,l,a_in,a_o);
    B_esl = MagneticField_ESL(x_d,y_d,z_d,l,Je_dw,ae);
    Bx_ml(j) = B_ml(1);
    Bx_esl(j) = B_esl(1);
    Bz_ml(j) = B_ml(3);
    Bz_esl(j) = B_esl(3);
    
    Bx_sum = 0;
    By_sum = 0;
    Bz_sum = 0;

    for i = 1:np
        
        R = (sqrt((x_d-p(1,i)).^2 + (y_d-p(2,i)).^2 + (z_d-p(3,i)).^2)).^3;

        Bx = 1*m(i)/(4*pi)*(x_d-p(1,i))./R;
        By = 1*m(i)/(4*pi)*(y_d-p(2,i))./R;
        Bz = 1*m(i)/(4*pi)*(z_d-p(3,i))./R;

        Bx_sum = Bx_sum + Bx;
        By_sum = By_sum + By;
        Bz_sum = Bz_sum + Bz;
    end
    Bx_dmp(j) = nu0*Bx_sum;
    Bz_dmp(j) = nu0*Bz_sum;
    
end
plot(r,-J*Bx_ml)
% plot(r/a_o,-J*Bz_ml,'r',r/a_o,-Bz_esl,'b',r/a_o,J*Bz_dmp,'k')

% plot(r/a_o,-J*Bx_ml,'r','LineWidth',3)
% hold on
% plot(r/a_o,-Bx_esl,'b--','LineWidth',2)
% plot(r/a_o,J*Bx_dmp,'k','LineWidth',2)

% plot(r/a_o,-J*Bz_ml,'r','LineWidth',3)
% hold on
% plot(r/a_o,-Bz_esl,'b--','LineWidth',2)
% plot(r/a_o,J*Bz_dmp,'k','LineWidth',2)

% plot(2*r/l,-J*Bz_ml,'r','LineWidth',3)
% hold on
% plot(2*r/l,-Bz_esl,'b--','LineWidth',2)
% plot(2*r/l,J*Bz_dmp,'k-','LineWidth',2)
% -J*Bz_ml(1)
% -Bz_esl(1)
%%Normalized
% plot(r/a_o,Bz_ml/Bz_ml(1),'r','LineWidth',3)
% hold on
% plot(r/a_o,Bz_esl/Bz_esl(1),'b--','LineWidth',2)
% plot(r/a_o,Bz_dmp/Bz_dmp(1),'k','LineWidth',2)
xlabel('2z/l')
% xlabel('r/a_o')
ylabel('B_z (T)')
legend('ML','ESL','DMP_M_L')
% axis([1,6,0,0.03])
    